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Abstract.  We  study  the  nonlinear  dynamics  of  a Penning  trap  plasma,  including  the  effect  of 
the  finite  length  and  end  curvature  of  the  plasma  column.  A new  cylindrical  PIC  code,  called 
KANDINSKY,  has  been  implemented  by  using  a new  interpolation  scheme.  The  principal  idea  is  to 
calculate  the  volume  of  each  cell  from  a particle  volume,  in  the  same  manner  as  it  is  done  for  the 
cell  charge.  With  this  new  method,  the  density  is  conserved  along  streamlines  and  artificial  sources 
of  compressibility  are  avoided.  The  code  has  been  validated  with  a reference  Eulerian  fluid  code. 
We  compare  the  dynamics  of  three  different  models:  a model  with  compression  effects,  the  standard 
Euler  model  and  a geophysical  fluid  dynamics  model.  The  results  of  our  investigation  prove  that 
Penning  traps  can  really  be  used  to  simulate  geophysical  fluids. 


INTRODUCTION 

Malmberg-Penning  traps  are  often  characterized  as  a tool  to  simulate,  in  a plasma,  the  2D 
dynamics  described  by  the  Euler  equations.  However,  it  is  well  known  that  the  classical 
2D  drift-Poisson  model  (mathematically  equivalent  to  the  2D  Euler  model  [ 1 ])  is  unable 
to  capture  the  features  of  the  k = 1 diocotron  instability.  In  this  approximation,  there  are 
no  unstable  modes  [2]  and  the  continuum  spectrum  can  give  only  algebraic  growth  [3]; 
instead,  experiments  show  an  exponential  growth  of  the  mode  [4], 

Recently,  the  research  has  focused  on  trying  to  resolve  this  contradiction  between 
theory  and  experiments.  Finn  el  al.  [5]  were  the  first  to  explain  the  instability  in  terms 
of  effects  due  to  the  finite  length  and  end  curvature  of  the  plasma  column  (compression 
effects)  and  later  Coppa  et  al.  [6]  improved  this  fluid  model  by  modifying  the  expression 
for  the  velocity  field,  by  including  temperature  effects  and  by  using  the  exact  Green’s 
function  for  the  plasma  length.  Hilsabeck  and  O’Neil  [7]  pointed  out  that  kinetic  effects 
might  also  be  important  in  the  k = 1 diocotron  instability. 

Interestingly,  when  compression  effects  are  taken  into  account,  the  classical  analogy 
between  the  2D  drift-Poisson  model  and  the  Euler  model  for  an  incompressible  and 
inviscid  fluid  breaks  down.  Instead,  a new  analogy  between  non-neutral  plasmas  and 
geophysical  fluid  dynamics  (GFD)  arises  [5,  8],  in  which  the  plasma  length  plays  the 
role  of  the  inverse  of  the  fluid  depth.  This  analogy  is  very  important  from  the  point  of 
view  of  performing  geophysical  experiments  in  Penning  traps. 

In  the  present  contribution,  we  focus  on  a simplified  version  of  the  model  presented  in 
Ref.  [6],  obtained  by  neglecting  the  temporal  variation  of  the  plasma  length.  As  shown 
in  Ref.  [6],  this  model  is  still  k = 1 unstable.  The  model  is  compared  with  the  Euler 
model  and  with  the  y-plane  approximation  of  the  geophysical  fluid  dynamics. 

We  present  simulations  of  the  nonlinear  dynamics  obtained  with  the  cylindrical  PIC 
Code  KANDINSKY,  showing  that  GFD  can  be  actually  simulated  in  Penning  traps. 
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PHYSICAL  MODEL 


A cylindrical  Malmberg-Penning  trap  confining  an  electron  plasma  is  considered  in  the 
following.  The  trap  consists  of  three  electrodes.  The  central  cylinder,  which  extends 
between  z = — Lc/2  and  z = Lc/ 2,  is  grounded,  while  the  two  cylindrical  end  caps  are  at 
a negative  potential,  -V.  An  uniform  magnetic  field  provides  radial  confinement,  while 
the  potential  difference  between  the  central  electrode  and  the  end  caps  provides  axial 
confinement. 

In  typical  experiments,  the  electron  cyclotron  radius  is  much  smaller  than  the  typical 
size  of  the  device,  allowing  the  guiding  center  approximation.  Furthermore,  the  depen- 
dence of  the  plasma  properties  along  the  axial  direction  can  be  simplified  considerably 
thanks  to  the  rapid  bouncing  motion  between  the  end  caps.  However,  a correct  descrip- 
tion of  the  diocotron  instability,  especially  for  the  k = 1 mode,  requires  one  to  include 
the  effect  of  finite  length  of  the  plasma  in  the  axial  direction  and  end  curvature. 

Following  the  approach  in  Refs.  [5,  6],  the  particles  are  described  by  strings  of 
variable  length,  which  change  their  axial  length  and,  consequently,  their  density,  if 
they  move  radially  due  to  the  E x B drift.  The  effective  axial  length  is  calculated  self- 
consistently  by  solving  the  Poisson  equation  in  the  trap  and  assuming  thermodynamic 
equilibrium  along  the  axial  direction  [9],  In  normalized  units,  the  model  consists  of  the 
following  system  of  equations  [6]: 


< 


<t>eff  = <t>c  + alog 


L 

Tr— 0 


(1) 


where  a(r,9,f) 


o is  the  density  integrated  along  axial  direction,  Vx  (r,0,?)  = 

£o<tVo(0) 


is  the  plasma  velocity  in  the  plane  (r,8),  L(r,  0,r)  = ^ is  the  effective  length 
of  the  plasma,  T'(r,0,z,r)  = is  the  correction  to  the  potential  in  a trap 

<t>rO(0) 

of  finite  length  with  respect  to  the  value  in  a trap  of  infinite  length,  <)>eff(ri  9,  t)  is  the 

effective  potential,  a = is  the  electron  temperature,  t = ^^?is  the  time  and  Rw  is 

the  wall  radius.  Hatted  quantities  are  physical  and  have  dimensions;  the  corresponding 
normalized  quantities  share  the  same  symbol  but  are  not  hatted.  Subscripts  c and  0 label, 
respectively,  quantities  evaluated  in  the  central  plane  and  unperturbed  quantities.  For 
further  details  on  the  model  and  its  derivation  we  refer  to  Refs.  [6]  and  [10]. 

The  model  can  be  regarded  as  quasi-2D  since  only  the  potential  in  the  trap  depends 
on  the  axial  coordinate  and  can  be  reduced  to  a fully  2D  model  by  neglecting  the  time 
variation  of  the  effective  plasma  length  (in  this  case  one  does  not  need  to  solve  the  3D 
Poisson  equation).  In  this  approximation  the  new  model  becomes: 
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< 


(2) 


|(/J,Zo)  = -V1-V1KZo) 
v±  = e-  X V_L  (J)cff 
Vl  (J)c  = «c 

/q 

4>cff  = <t>c  + alOg- 

-MJ,r=0 


where  we  used  the  fact  that  a = ncLo  and  Lq  ( r ) is  obtained  from  the  initial  equilibrium 
solution. 

Compression  effects  due  to  the  finite  length  of  the  plasma  are  included  in  the  first  and 
second  equation  of  system  (2).  In  the  limit  of  a = 0 and  To/To, r=o  = 1 - sr  the  model 
described  by  Eqs.  (2)  reduces  to  the  model  proposed  by  Finn  et  al.  (Eq.  (38),  Ref.  [5]a). 

In  the  following  we  will  focus  our  attention  on  model  (2),  since  it  is  simplified  but 
still  able  to  capture  the  features  of  the  k = I diocotron  instability  (as  shown  in  Ref.  [6]). 
The  spectrum  of  eigenvalues  of  model  (2)  is  investigated  in  Ref.  [11], 

Our  aim  is  to  compare  the  dynamics  of  model  (2)  with  the  standard  Euler  model: 


i 


dnc 


= -' Vi-Vja- 
— e-  xVj  <)v 


Vj  <tv  = nc 


and  with  the  geophysical  fluid  dynamics  (GFD)  model  [12]: 


(3) 


d_ 

dt  V h 
Vi  = e;  x V±  *F 


-V.  -V, 


(4) 


where  E,  is  the  vonicity,  Q is  the  rotation  frequency,  h is  the  height  of  the  free  surface  of 
the  fluid  relative  to  the  bottom  topography,'}'  is  the  strcamfunction  and  (£  + 2£2)//i  is  the 
potential  vorticity.  Model  (4)  is  obtained  from  the  shallow-water  equations  in  the  limit  of 
small  Rossby  number,  namely  in  the  quasi-geostrophic  approximation.  As  pointed  out 
in  Ref.  [5],  with  = 0 the  effective  length  of  the  plasma  To  plays  in  model  (2)  the  same 
role  played  by  l//i  in  GFD.  It  is  this  analogy  which  leads  to  the  possibility  of  simulating 
experimentally  geophysical  fluid  dynamics  in  a Penning  trap. 

In  all  the  models,  a quantity  is  advected  by  the  velocity  field:  the  density  (vorticity) 
tic  for  Euler,  the  line  integrated  density  ncL o for  model  (2)  and  the  potential  vorticity 
for  GFD.  All  these  models  have  the  same  equation  for  the  velocity  field,  but  only  for  the 
Euler  case  is  the  velocity  field  directly  calculated  from  the  quantity  which  is  advected.  In 
model  (2)  the  velocity  field  is  calculated  from  the  solution  of  the  Poisson  equation  plus 
an  extra  term  which  retains  curvature  and  thermal  effects,  while  in  GFD  the  velocity 
field  is  determined  by  the  vorticity  £,. 


PIC  CODE 

The  PIC  code  KANDINSKY  has  been  developed  based  on  a standard  cylindrical  PIC 
code  previously  developed  [13]  by  the  authors  in  collaboration  with  Gianni  Coppa 
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and  Antonio  D’ Angola  at  the  Politecnico  di  Torino.  A cylindrical  PIC  code  has  many 
advantages:  very  low  viscosity;  boundary  conditions  can  be  imposed  exactly;  the  Fast 
Fourier  Transform  (FFT)  can  be  used  to  solve  the  Poisson  equation  resulting  in  an 
improvement  of  the  speed  of  the  code.  Nevertheless,  it  is  well  known  that  cylindrical 
PIC  codes  can  potentially  suffer  from  a number  of  problems.  Of  particular  importance 
to  the  present  study  are  the  noise  that  typically  arises  in  the  center  of  the  system,  the  non- 
conservation of  enstrophy,  and  the  growth  of  the  maximum  value  of  the  density  during 
simulations  (inconsistent  with  incompressible  flows).  In  the  following  subsections  we 
will  show  how  to  improve  the  cylindrical  PIC  scheme  and  to  solve  the  problem. 


The  problem 

In  order  to  understand  the  origin  of  the  problems  that  affect  cylindrical  PIC  codes,  we 
will  start  with  a simple  exercise.  Consider  the  ID  uniform  cartesian  grid  of  Fig.  la  and 
consider  one  particle  of  unit  charge  located  at  positions  A,  B or  C.  It  is  possible  to  cal- 
culate the  density  of  each  cell  with  the  volume-weighted  (Cloud-In-Cell)  interpolation 
scheme,  which  is  the  most  common  scheme  adopted  in  PIC  codes.  When  the  particle 
is  in  position  A,  it  contributes  only  to  the  first  cell  and  the  density  is  pi  = 1/Ax  (Ax  is 
the  grid  spacing);  when  the  particle  is  in  position  B,  it  is  spread  equally  on  the  two  cells 
and  pi  = p2  = l/(2Ax);  when  the  particle  is  in  position  C,  it  contributes  only  to  the 
second  cell  and  pi  — 1 /Ax.  During  the  motion  of  the  particle  from  A to  C,  the  density 
at  the  particle  position  fluctuates  between  l/(2Ax)  and  I /Ax  but  on  average  it  is  con- 
served over  long  Lagrangian  trajectories.  Furthermore,  this  fluctuation  is  reduced  when 
summed  over  the  number  of  particles  per  cell,  and  this  source  of  noise  is  manageable. 

We  now  repeat  the  same  exercise  for  the  uniform  cylindrical  grid  of  Fig.  lb.  In  general 
a particle  would  contribute  to  two  cells  in  0 and  two  cells  in  r direction,  since  the  grid  is 
2D.  In  this  specific  example,  we  place  A,  B and  C at  the  same  azimuthal  angle  A0/2,  so 
there  is  only  contribution  to  one  cell  in  0 and  the  problem  can  be  considered  ID.  When 
the  particle  is  in  position  A,  it  contributes  only  to  the  first  cell  and  pi  = 2/(Ar2A0)  (Ar 
and  A0  are,  respectively,  the  grid  spacing  in  r and  0 directions);  when  the  particle  is  in 
position  B,  its  charge  is  spread  equally  between  the  two  cells  but  the  different  volume  of 
the  two  cells  leads  to  different  densities,  pi  = l/(Ar2A0)  and  p2  = l/^A/^Ad).  When 
the  particle  is  in  position  C,  it  contributes  only  to  the  second  cell  and  P2  = 2/ (3Ar2AQ). 
Again  during  the  motion  from  A to  C the  density  is  transferred  from  the  first  cell  to  the 
second,  but  in  this  case  the  total  density  of  the  particle  does  not  remain  constant,  even  on 
average.  In  particular,  the  density  decreases  as  the  particle  moves  away  from  the  center  of 
the  system  and  this  means  that  the  CIC  interpolation  scheme,  applied  to  a cylindrical  grid 
with  uniform  spacings  Ar  and  A0,  intrinsically  violate  the  incompressibility  condition. 
This  effect  is  partially  compensated  for  on  average  because  more  particles  are  in  the  cell 
of  volume  rArA0  for  large  r. 

This  simple  exercise  points  out  that,  when  dealing  with  non-uniform  cell  volumes, 
standard  interpolation  schemes  produce  fictitious  compressibility. 

We  emphasize  that  this  problem  affects  not  only  the  cylindrical  geometry  but  it  occurs 
whenever  a grid  of  non-uniform  volume  is  used. 
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FIGURE  1.  a)  Uniform  cartesian  grid  (left);  b)  Uniform  cylindrical  grid  (right). 


The  solution 

According  to  [14],  we  define  natural  coordinates,  to  = (C,r|)  by  mapping  each  cell  of 
the  physical  grid  onto  a unit  square  cell  in  the  space  of  natural  coordinates.  Therefore,  in 
the  space  of  natural  coordinates  the  grid  is  uniform  and  cartesian.  We  define  a shape 
function  .?(£>  - <s)p)  in  order  to  interpolate  from  the  particles  to  the  grid  and  vice- 
versa,  normalized  to  / f sdCjd Tj  = 1 in  the  space  of  natural  coordinates.  Then,  the  shape 
function  in  the  physical  space  is  S(x( to)  - xP(&p))  = s(oi  - &p).  When  the  grid  and  the 
particles  move  with  the  fluid  velocity  and  the  interpolation  scheme  is  consistent  with 
the  choice  of  the  grid,  it  can  be  shown  that  the  natural  coordinates  are  constants  of  the 
motion  [14];  in  other  words,  s(&  - o'),,)  is  constant  during  the  Lagrangian  phase  in  the 
space  of  natural  coordinates. 

In  order  to  solve  the  problem  that  we  described  in  the  previous  paragraph,  we  in- 
troduce the  volume  associated  with  each  particle,  Vp,  in  the  same  manner  as  it  is  done 
lor  the  particle  charge  qp,  and  to  calculate  the  volume  of  each  cell  we  use  the  same 
interpolation  scheme  adopted  for  the  charge: 

V(<3)  = - (7V)>  (5) 

p 

where  w is  the  assignment  function  associated  with  the  shape  function  ,v  [15].  Then  the 
density  is  given  by: 


p(d>)  - 


lpq,M&-Cap) 

lpVpw(C)  - (hp)' 


(6) 


This  definition  gives  conservation  of  p along  the  Lagrangian  trajectories  of  each  particle, 
since  qp  and  Vp  are  constant  (determined  at  the  beginning  of  each  simulation)  and  each 
particle  is  advected  with  constant  w in  the  space  of  natural  coordinates. 


PIC  summary 


The  KANDINSKY  code  consists  of  four  parts:  the  loading  of  particles,  the  interpo- 
lation scheme,  the  Poisson  solver  and  the  particle  mover.  The  implementation  of  each 
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part  of  the  code  has  been  realized  taking  into  account  the  characteristics  of  the  physical 
system  to  be  simulated. 

• Loading  particles 

In  order  to  provide  a precise,  noiseless  representation  of  the  initial  electron  density 
distribution,  computational  particles  are  allowed  to  have  different  charge.  This  can  be 
obtained  by  using  the  Mass  Matrix  formulation  [16],  according  to  which  the  particle 
properties  are  calculated  at  the  beginning  of  each  simulation. 

• Interpolation  particle-grid 

The  density  on  the  grid  is  calculated  from  the  particles  through  the  new  interpola- 
tion scheme  defined  in  Eq.  (6).  The  classic  Cloud-In-Cell  (CIC)  method  [15]  is  used 
for  the  interpolation  functions.  Slightly  different  expressions  are  used  when  a particle 
approaches  the  center  or  the  boundary  of  the  domain. 

• Poisson  solver 

The  Poisson  equation  can  be  solved  with  an  efficient  algorithm  based  upon  the  Fast 
Fourier  Transform  in  0,  assuming  a piecewise  constant  density  distribution  in  r [13].  A 
set  of  decoupled  differential  equations  is  obtained  for  the  Fourier  components  <)>*(/-): 


I ±(r  4A  + 1 

r dr  \ dr  j r2 


2 

$k=nk(r). 


(7) 


Eqs.  (7)  are  solved  using  a finite  difference  scheme.  Then  the  velocity  field  can  be 
calculated  by  using  the  central  difference  discretization. 

• Interpolation  grid-particle 

Given  the  velocity  field  on  the  grid,  the  velocity  on  each  particle  can  be  calculated 
using  the  same  CIC  interpolation  scheme: 

VP  = ^vcw(d)c-tip).  (8) 

C 

• Particle  mover 

Particles  are  moved  solving  the  equations  of  motion  drp/dt  = vp  discretized  with  the 
fourth-order  Runge-Kutta  method. 

Figures  2 and  3 compare  two  simulations  that  only  differ  in  using  (6)  (in  the  following 
we  will  refer  to  this  code  as  an  Advanced  PIC  code  or  KANDINSKY  code)  or  the  usual 
P = YjPcipwl (rArAS)  (Standard  PIC  code)  for  interpolating  the  density  on  the  grid.  This 
run  has  been  done  for  an  initial  hollow  density  profile  with  a ft  = 2 perturbation.  It  can 
be  noticed  that  very  early  some  differences  arise,  in  particular  the  Standard  PIC  code 
develops  some  high  density  vortices  (the  maximum  value  of  the  density  almost  doubles 
from  the  first  to  the  last  picture  of  Fig.  2)  while  this  does  not  happen  for  the  Advanced 
PIC  code  (the  maximum  of  the  density  remains  constant).  In  the  Standard  PIC  code, 
we  add  an  artificial  compressibility  to  the  system  which  is  manifested  by  the  localized 
vortices  with  increasing  density. 

Furthermore,  in  the  Standard  PIC  code  enstrophy  grows  because  of  the  contribution 
of  high  density  vortices,  while  we  would  expect  enstrophy  to  be  constant  (as  it  is  theoret- 
ically for  the  Euler  model)  or  to  decrease  slowly  (because  the  density  develops  smaller 
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FIGURE  2.  Standard  PIC  Code  for  a k = 2 instability  in  a hollow  profile. 
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FIGURE  3.  KANDINSKY  Code  for  a Jt  = 2 instability  in  a hollow  profile. 


and  smaller  length  scales  and  at  some  point  it  is  dissipated  by  numerical  viscosity). 
Indeed,  for  the  Advanced  PIC  code,  enstrophy  is  a slowly  decreasing  function  of  time. 


VALIDATION  TESTS 

In  this  section  we  use  the  Euler  model,  which  can  be  obtained  from  system  (2)  by  taking 
Lo  (r)  =const,  in  order  to  compare  our  results  with  analytical  results  and  with  a reference 
Eulerian  fluid  code.  In  all  of  our  simulations  wc  use  an  initial  equilibrium  density  profile 
of  the  form 


r 

/ r \ 2" 

2 

/ r \ 2 

nco(r)  = < 

1 /»0  (0) 
| 

,+,"+2»U). 

l o 

0 < r < rp 
r>rp , 


(9) 


which  depends  on  three  parameters:  «o(0),  rp  and  The  last  one  is  particularly  im- 
portant: for  n > 0 the  profile  is  hollow  (possible  diocotron  instability)  where  for  g < 0 
the  profile  is  monotonic  (always  stable).  The  simulations  below  are  done  with  /t  = 10, 
rp  = 0.59  and  no  (0)  is  specified  to  give  J'q"  nta(r)2nrdr  = 1 . The  normalized  tempera- 
ture a is  chosen  equal  to  0.  According  to  (6]  the  equilibrium  profile  of  a is  obtained  by 
Go  = HcqLo-  The  initial  charge  density  is  perturbed  by  modifying  the  particle  charge  as 
q'p  = qp  [1  + eco.v(£0)],  where  k is  the  mode  number  and  e is  the  perturbation  amplitude. 
The  system  is  discretized  with  a Nr  x No  grid  and  using  Np  particles.  In  the  cases  below 
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Vorticity  Perturbation  at  Jim©  12.5 


FIGURE  4.  Vorticity  (density)  perturbation  from  a)  Eulerian  Fluid  Code:  linear  phase  (left);  b) 
KANDINSKY  Code:  linear  phase  (right). 


we  use  Nr  = 150,  Nq  = 128,  Np  = 558208,  the  average  number  of  particles  per  cell  is  49 
and  the  time  step  is  dt  = 0. 1. 

• Test  of  accuracy 

Three  important  conserved  quantities  are  canonical  angular  momentum,  propor- 
tional to  P = f f i^nc-rdrdQ,  electrostatic  energy,  proportional  to  H — J J nc(\>crdrdQ, 
and  enstrophy,  proportional  to  Z — J f nj.rdrd0.  (There  are  additional  invariants 
f 1 f(nc)rdrd 0,  with  / any  function.) 

The  first  two  quantities  are  conserved  well  by  the  code  (in  a run  with  k = 2,  e = 0.01 
and  tmax  = 600,  not  shown):  8P/P  = 0.6%  and  8H/H  = 0.3%.  The  enstrophy  is  weakly 
conserved  due  to  filamentation,  but  still  the  code  preserved  it  well:  8 Z/Z  = 5%. 

• Comparison  with  a Reference  Eulerian  Fluid  Code 

The  code  has  been  compared  with  a reference  Eulerian  fluid  code  in  both  the  linear 
and  nonlinear  phase.  Figs.  4a  and  4b  show  the  evolution  of  the  perturbation  of  the  density 
at  t = 12.5.  The  initial  perturbation  has  been  chosen  very  small  (e  = 0.001)  in  order  to 
remain  in  the  linear  phase.  The  agreement  between  the  two  codes  is  remarkable  and  one 
can  notice  that  the  fluid  code  has  more  numerical  viscosity  than  the  PIC  code. 

Figs.  5a  and  5b  show  the  evolution  of  the  density  after  a large  initial  perturbation 
(e  = 0.25).  The  two  codes  agree  very  well,  even  in  this  nonlinear  case.  It  should  be 
noticed  that  the  Eulerian  code  cannot  follow  with  enough  accuracy  the  region  with  steep 
gradients  outside  the  core  region. 


SIMULATIONS 

In  this  section,  we  analyse  the  dynamics  of  models  (2),  (3)  and  (4).  We  compare 
different  runs  with  k — 1 and  £ = 0.001.  Model  (2)  can  simulate  GFD  by  using  a 
suitable  equilibrium  length  profile,  as  it  will  be  explained  in  the  following.  The  shallow- 
water  equation  for  the  potential  vorticity  can  be  simplified  assuming  small  topography 
variations  or  small  variations  of  the  Coriolis  parameter.  In  the  last  case,  (£  + 2£2)//j  can 
be  approximated  as  ^ + (3r  — yr2.  Fory  = 0 we  have  the  (3-plane  approximation,  while 
for  (3  = 0 we  have  the  y-plane  approximation  [12],  The  analogy  between  non-neutral 
plasmas  and  GFD  is  established  using  the  linear  theory,  where  3 — 2y r can  be  identified 
with  nco(r)£^(r) / Lo(r)  [5,  8],  Moreover,  due  to  the  fact  that  £^(0)  = 0,  it  follows  that 
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Vorboty  <1  l«n«  6 20 


FIGURE  S.  Vorticity  (density)  from  a)  Hitlerian  Fluid  Code:  nonlinear  phase  (left);  b)  KANDINSKY 
Code:  nonlinear  phase  (right). 


Vorway  mm  wm  Vo«ik*,  mm  36  S3  turn  Vw*Xy  mm  57  S3  turr* 


FIGURE  6.  Evolution  of  the  density  ii,  according  to  the  Euler  model. 


P = 0 and  thus  we  are  in  the  limit  of  the  y-plane  approximation.  Therefore,  in  order  to 
simulate  the  y-plane  approximation  of  GFD,  we  use  model  (2)  with  a = 0 and  the  length 
profile  obtained  by  solving 

= -2yr,  (10) 

with  the  initial  density  profile  given  by  (9)  and  y = .v»o(0)  (this  is  the  value  of  y 
obtained  from  (10)  assuming  a parabolic  profile  Lq  — 1 - .vr2,  a constant  density  profile 
iicQ  — ;io(0)  and  in  the  limit  of  small  radius,  as  shown  in  Ref.  [8]).  Model  (4)  is  then  also 
simulated  in  the  limit  of  y-plane  approximation,  using  the  same  y = sno{0).  In  all  the 
simulations  below  we  use  the  curvature  s = 0.5. 

Figure  6 shows  the  evolution  of  the  density  for  the  Euler  model  after  an  initial 
k = 1 perturbation.  The  model  does  not  develop  any  k = 1 instability  and  a slowly 
growing  k = 2 mode  eventually  dominates.  Figure  7 shows  the  density  at  the  same  times 
according  to  model  (2).  This  model  exhibits  a strong  k = 1 diocotron  instability  which 
dominates  the  k = 2 mode  seen  in  Fig.  6.  The  low-density  core  region  is  pushed  off-axis 
and  the  high-density  region  fills  the  center  of  the  trap.  At  later  times  (not  shown),  the 
low-density  region  is  redistributed  by  fi lamentation  or  in  big  vortices.  The  same  picture 
can  be  drawn  from  Fig.  8,  which  shows  the  vorticity  (density)  evolution  of  the  GFD 
model.  The  density  develops  the  same  features  in  Figs.  7 and  8 with  small  differences  at 
later  times  due  to  nonlinearities. 

Clearly,  the  comparison  of  Figs.  6,  7 and  8 proves  that  a Penning  trap  can  simulate 
GFD  and  particularly  its  peculiar  features  not  present  in  the  classic  Euler  model. 
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Vorticrty  after  28  89  turns 


mode  k *»  1.  £■  0.001 


Vortidty  after  386  turns 


-0.5  0 0.5 

mode  k - 1.  £-  0.001 


Vorlioty  after  58.76  turns 


-0.5  0 0.5 

modek  ■ 1.  £=  0.001 


FIGURE  7.  Evolution  of  the  density  n,  according  to  model  (2). 


Vortidty  after  28.94  turns 


Vortidty  after  38.67  turns 


mode  k-  1,  £■  0 001 


Vortidty  after  58.89  turns 


-0  5 0 0.5 

modek  - 1.  £» 0.001 


FIGURE  8.  Evolution  of  the  vorticity  £,  according  to  the  GFD  model. 
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